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Abstract 


ALICE is one of four large experiments at the CERN Large Hadron Collider near Geneva, specially 
designed to study particle production in ultra-relativistic heavy-ion collisions. Located 52 meters 
underground with 28 meters of overburden rock, it has also been used to detect muons produced by 
cosmic ray interactions in the upper atmosphere. In this paper, we present the multiplicity distribu¬ 
tion of these atmospheric muons and its comparison with Monte Carlo simulations. This analysis 
exploits the large size and excellent tracking capability of the ALICE Time Projection Chamber. 
A special emphasis is given to the study of high multiplicity events containing more than 100 re¬ 
constructed muons and corresponding to a muon areal density p u > 5.9 m ' 2 . Similar events have 
been studied in previous underground experiments such as ALEPH and DELPHI at LEP. While these 
experiments were able to reproduce the measured muon multiplicity distribution with Monte Carlo 
simulations at low and intermediate multiplicities, their simulations failed to describe the frequency 
of the highest multiplicity events. In this work we show that the high multiplicity events observed in 
ALICE stem from primary cosmic rays with energies above 10 16 eV and that the frequency of these 
events can be successfully described by assuming a heavy mass composition of primary cosmic rays 
in this energy range. The development of the resulting air showers was simulated using the latest 
version of QGSJET to model hadronic interactions. This observation places significant constraints 
on alternative, more exotic, production mechanisms for these events. 
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1 Introduction 

ALICE (A Large Ion Collider Experiment) |l|] designed to study Quark-Gluon Plasma (QGP) formation 
in ultra-relativistic heavy-ion collisions at the CERN Large Hadron Collider (LHC), has also been used to 
perform studies that are of relevance to astro-particle physics. The use of high-energy physics detectors 
for cosmic ray physics was pioneered by ALEPH H, DELPHI il and L3 [4] during the Large Electron- 
Positron (LEP) collider era at CERN. An extension of these earlier studies is now possible at the LHC, 
where experiments can operate under stable conditions for many years. ALICE undertook a programme 
of cosmic ray data taking between 2010 and 2013 during pauses in collider operations when there was 
no beam circulating in the LHC. 


Cosmic ray muons are created in Extensive Air Showers (EAS) following the interaction of cosmic ray 
primaries (protons and heavier nuclei) with nuclei in the upper atmosphere. Primary cosmic rays span 
a broad energy range, starting at approximately 10 9 eV and extending to more than 10 2() eV. In this 
study, we find that events containing more than four reconstructed muons in the ALICE Time Projection 
Chamber (TPC), which we refer to as multi-muon events, stem from primaries with energy E > 10 14 eV. 
The detection of EAS originating from interactions above this energy, in particular around the energy of 
the knee in the primary spectrum ( E^ ~ 3 x 10 15 eV), has been performed by several large-area arrays at 
ground level (e.g. [5M3]), while deep underground detectors (e.g. |[8- J_C ]) have studied the high energy 
muonic component of EAS. The main aims of these experiments were to explore the mass composition 
and energy spectrum of primary cosmic rays. 

The muon multiplicity distribution (MMD) was measured at LEP with the ALEPH detector & This 
study concluded that the bulk of the data can be successfully described using standard hadronic produc¬ 
tion mechanisms, but that the highest multiplicity events, containing around 75-150 muons, occur with 
a frequency which is almost an order of magnitude above expectation, even when assuming that the pri¬ 
mary cosmic rays are purely composed of iron nuclei. A similar study was carried out with the DELPHI 
detector, which also found that Monte Carlo simulations were unable to account for the abundance of 
high muon multiplicity events [12]. Several proposals have been put forward in the scientific literature to 
explain this discrepancy. Some authors suggest that hypothetical strangelets form a small percentage of 
very energetic cosmic rays [13]], while others have tried to explain the excess of high muon multiplicity 
events by the creation of the QGP in interactions involving high mass primary cosmic rays (iron nuclei) 
with nuclei in the atmosphere [14]. 


In this paper, we exploit the large size and excellent tracking capability of the ALICE TPC [15Q to 
study the muonic component of EAS. We describe the analysis of the muon multiplicity distribution with 
particular emphasis on high muon multiplicity events containing more than 100 muons in a single event 
and corresponding to an areal density p u > 5.9 m 2 . We employ a description of the shower based upon 
the latest version of QGSJET [ 16, 17], a hadronic interaction model commonly used in EAS simulations. 


Details of the environment of ALICE and the detectors used for this analysis are described in the follow¬ 
ing section, while the selection of the data and the algorithm adopted to reconstruct atmospheric muons 
are discussed in Section 3. The muon multiplicity distribution and the study of high muon multiplicity 
events are described in Section 4. The results are presented in Section 5 and, finally, in Section 6 we 
make some concluding remarks. 


2 The ALICE experiment 

ALICE is located at Point 2 of the LHC accelerator complex, approximately 450 nr above sea level in a 
cavern 52 nr underground with 28 nr of overburden rock. The rock absorbs all of the electromagnetic and 
hadronic components of the observed EAS, so that only muons with an energy E, at the surface, larger 
than 16 GeV reach the detectors [180. The geometry of ALICE is typical of a collider experiment. A 
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large solenoidal magnet forms a central barrel that houses several detectors, including a large, cylindrical 
TPC. Outside the solenoid, and on one end, there is a single-arm, forward spectrometer, which was not 
used in this analysis. A complete description of the apparatus is given in 0 ]]. 

The ALICE TPC is the largest detector of its type ever built. It was used to reconstruct the trajectory 
of cosmic ray muons passing through the active volume of the detector, which comprises a cylindrical 
gas volume divided into two halves by a central membrane. The TPC has an inner radius of 80 cm, an 
outer radius of 280 cm and a total length of 500 cm along the LHC beam direction. At each end of 
the cylindrical volume there are multi-wire proportional chambers with pad readout. For the purpose 
of detecting cosmic ray muons, the total area of the detector due to its horizontal cylindrical geometry 
is approximately 26 nr 2 . However, after placing a cut on the minimum length required to reconstruct a 
cosmic ray track the maximum effective area reduces to approximately 17 nr 2 . The apparent area of the 
detector also varies with the zenith angle of the incident muons. Track selection is discussed in more 
detail in Section [3] An example of a single atmospheric muon event is shown in Fig. [T] 



Fig. 1: A single atmospheric muon event. The thin outer cylinder is the Time Of Flight detector (1). The large 
inner cylinder is the Time Projection Chamber (2) and the smaller cylinder at the centre is the silicon Inner Tracking 
System (3). Muons are reconstructed as two TPC tracks, one in the upper half of the detector (up track) and the 
other in the lower half (down track), which are then joined to create a single muon track. 


Three detector subsystems were used to provide dedicated triggers for this study: ACORDE (Alice 
COsnric Ray DEtector) [19], SPD (Silicon Pixel Detector) [20] and TOF (Time Of Flight detector) |21], 


ACORDE is an array of 60 scintillator modules located on the three upper faces of the octagonal yoke of 
the solenoid, covering 10% of its surface area. A trigger was formed by the coincidence of signals in two 
different modules (a two-fold coincidence), although the trigger can also be configured to select events 
when a single module fires or when more than two modules fire. 
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The SPD is part of the Inner Tracking System located inside the inner field cage of the TPC. It is com¬ 
posed of two layers of silicon pixel modules located at a distance of 39 mm and 76 mm from the LHC 
beam axis, respectively. The layers have an active length of 28.3 cm, centred upon the nominal inter¬ 
action point of the LHC beams. The SPD was incorporated into the trigger by requiring a coincidence 
between signals in the top and bottom halves of the outermost layer. 

The TOF is a cylindrical array of multi-gap resistive-plate chambers that completely surrounds the outer 
radius of the TPC. The TOF trigger requires a signal in a pad corresponding to a cluster of readout 
channels covering an area of 500 cm 2 in the upper part of the detector and another signal in a pad in the 
opposite lower part forming a back-to-back coincidence with respect to the central axis of the detector. 
The configuration of the pads involved in the trigger can be changed via software. In some periods of 
data taking, this flexibility has been exploited to require a signal in an upper pad and in the opposing pad 
plus the two adjacent pads forming a back-to-back ±1 coincidence. 

Cosmic ray data were acquired with a combination (logical OR) of at least two out of the three trigger 
conditions (ACORDE, SPD and TOF) depending on the run period. The trigger efficiency was studied 
with a detailed Monte Carlo simulation, which is discussed in Section |H Most events were classified as 
either single muon events or nrulti-nruon events, with a small percentage of “interaction” events where 
very energetic muons have interacted with the iron yoke of the magnet producing a shower of particles 
that pass through the TPC. 


3 Event reconstruction and data selection 


The TPC tracking algorithm [22] was designed to reconstruct tracks produced in the interaction region of 
the two LHC beams. It finds tracks by working inwards from the outer radius of the detector where, dur¬ 
ing collider operation, the Lack density is lowest. The present analysis used the same tracking algorithm 
but removed any requirement that tracks should pass through a central interaction point. However, the 
tracking algorithm has not been optimised for very inclined (quasi horizontal) tracks. Therefore, to avoid 
reconstruction inaccuracies associated with the most inclined showers, we restricted the zenith angle of 
all events to the range 0° < 6 < 50°. 


As a consequence of reconstructing Lacks from the outer radius of the TPC inwards, cosmic ray muons 
are typically reconstructed as two separate tracks in the upper and lower halves of the TPC as shown in 
Fig. |T] We refer to these tracks as up and down tracks. Following this first pass of the reconstruction a new 
algorithm was applied to match each up track with its corresponding down track to reconstruct the full 
trajectory of the muons and to eliminate double counting. Starting with single muon events (producing 
two TPC tracks), where the matching of tracks is straightforward, the reconstruction has been tuned to 
handle events containing hundreds of muons. High multiplicity Monte Carlo events have been used to 
optimise the matching performance. 


Each TPC Lack can be reconstructed with up to 159 individual space points. In order to maximise the 
detector acceptance for this analysis, tracks were required to have a minimum of 50 space points and, in 
events where the magnetic field was on, a momentum greater than 0.5 GeV/c to eliminate all possible 
background from electrons and positrons. In multi-muon events, accepted tracks were required to be 
approximately parallel since atmospheric muons coming from the same EAS arrive almost parallel at 
ground level. The parallelism cut involves forming the scalar product of the direction of the analysed 
track t a with a reference track t r , requiring that t a -t r = cos(A'P) > 0.990 to accept the analysed track. 
The reference track was chosen to give the largest number of tracks satisfying the parallelism cut. This 
requirement introduces an additional momentum cut due to the bending of muon tracks in the magnetic 
field. The momentum cut is a function of the azimuth angle of the muon track and varies between 1 
and 2 GeV/c. Finally, each up track was matched to the nearest down track if the distance of closest 
approach between them at the horizontal mid plane of the TPC was d xz < 3 cm. This value was chosen 
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to be large enough to maximise the matching efficiency in high multiplicity Monte Carlo events, while 
keeping combinatorial background to a minimum. 

A muon reconstructed with two TPC tracks (up and down) is called a “matched muon”. When a TPC 
track fulfils all the criteria to be a muon track: number of space points, momentum and parallelism, but 
does not have a corresponding track within d xz < 3 cm in the opposite side of the TPC, this track is still 
accepted as a muon candidate but flagged as a “single-track muon”. Most single-track muons are found 
to cross the TPC near its ends where part of the muon trajectory falls outside the detector. 

To quantify the performance of the tracking and matching algorithms, we studied the multiplicity depen¬ 
dence of the reconstruction efficiency using Monte Carlo simulated events. We generated 1000 events 
for 20 discrete values of the muon multiplicity, varying between 1 and 300, which were then recon¬ 
structed using the same algorithms applied to real events. In each event, muons were generated parallel 
to each other like in EAS and cross the whole TPC volume. Fig. [2] shows the mean values (MEAN) and 
root-mean-square (RMS) of the relative difference between the number of generated and reconstructed 
muons, 

(# generated muons — # reconstructed muons) / (# generated muons), (1) 

as a function of the number of generated muons. The root-mean-square represents the resolution on the 
number of reconstructed muons and is typically less than 4%, while for the highest multiplicities it is 
around 2%. The mean value is less than 1% up to /V tl rj 50, increasing to 5% at high muon multiplicities 
(Nfi ~ 300). 
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Fig. 2: Root-mean-square and mean values of the relative difference between the number of generated and recon¬ 
structed muons for events simulated with different muon multiplicities. 

To illustrate the similarity of the data and the Monte Carlo simulation, Fig.[3]shows the ratio of the num¬ 
ber of muons reconstructed as single tracks (either up or down tracks) to the total number of reconstructed 
muons (both single and matched tracks) for different multiplicities. The ratio obtained from the data is 
compared with the ratios obtained from simulated samples of pure proton primary cosmic rays and pure 
iron primaries. Over the range of intermediate muon multiplicities shown, the ratio varies between 0.2 
and 0.4 with good agreement between data and simulations. There is no significant difference between 
the simulated proton and non samples. 
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Fig. 3: The ratio of muons reconstructed as single tracks to the total number of reconstructed muons (both single 
and matched tracks) in the data and simulations with proton and iron primaries. 

Data were recorded between 2010 and 2013 during pauses in collider operations when no beam was 
circulating in the LHC. The total accumulated run time amounted to 30.8 days, resulting in approximately 
22.6 million events with at least one reconstructed muon (single-track or matched) in the TPC. Only 
multi-muon events are discussed further in this paper. We define multi-muon events as those events 
with more than four reconstructed muons in the TPC (N^ > 4). In total, we collected a sample of 7487 
multi-muon events. 

4 Analysis of the data and simulation 

To obtain the MMD we have corrected the measured distribution for the efficiency of the trigger. The 
correction was calculated from a Monte Carlo simulation that is described later in this section. Given 
the complementary coverage of the TOF barrel to the TPC, the TOF trigger was mainly responsible for 
selecting events in the low-to-intermediate range of muon multiplicities (7 < < 70). The efficiency 

of the TOF trigger as a function of the muon multiplicity is shown in Fig. [4] The efficiency is lower at low 
muon multiplicity due to the back-to-back coincidence requirement of the TOF trigger. The efficiency 
of the ACORDE trigger has a similar, increasing trend with the muon multiplicity. The multiplicities 
at which the two triggers reach full (100%) efficiency are /V u > 10 (TOF) and V tl > 15 (ACORDE). 
Given the much smaller area of the SPD in comparison with the TPC, the efficiency of the SPD trigger 
is significantly lower than both ACORDE and TOF. It makes only a minor contribution to the MMD in 
the low-to-intermediate range of muon multiplicities. 

The MMD obtained from the whole data sample and corrected for trigger efficiency is shown in Fig. [5] 
Values for the systematic uncertainty in the number of events as a function of multiplicity have been 
estimated by varying the parameters of the track reconstruction and matching algorithms. We find a 
smooth distribution up to a muon multiplicity of around 70 and then 5 events with a muon multiplicity 
greater than 100. We define the events with /V u > 100 high muon multiplicity (HMM) events. Given 
the nature and topology of high multiplicity events, all trigger conditions contributed to this sample with 
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Fig. 4: TOF trigger efficiency as a function of muon multiplicity. 

close to 100% efficiency. The aim of the following analysis is to model the MMD at low-to-intermediate 
multiplicities and to explore the origin of the HMM events. 



Fig. 5 : Muon multiplicity distribution of the whole sample of data (2010-2013) corresponding to 30.8 days of data 
taking. 


The difficulty in describing EAS, and consequently the number of muons reaching ground level, mainly 
arises from uncertainties in the properties of multi-particle production in hadron-air interactions. These 
interactions are often described phenomenologically within Monte Carlo event generators. Model param¬ 
eters, such as total and inelastic hadron-proton cross sections, inelastic scattering slopes and diffractive 
structure functions, are constrained by measurements obtained from accelerator experiments. 


In this analysis we have adopted the CORSIKA [23] event generator incorporating QGSJET [160 for 
the hadronic interaction model to simulate the generation and development of EAS. CORSIKA version 
6990 incorporating QGSJET 11-03 has been used to study the MMD distribution and HMM events; 
CORSIKA version 7350 incorporating QGSJET 11-04 has been used to check and confirm the results 
for HMM events. The significant differences between the two versions of QGSJET are the inclusion of 
Ponreron loops in the formalism of QGSJET 11-04 and a retuning of the model parameters using early 
LHC data for the first time [24]. Most relevant to the present study is that pion exchange is assumed to 
dominate forward neutral hadron production in the QGSJET 11-04, which has been shown to enhance the 
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production of p mesons resulting in an enhancement of the muon content of EAS by about 20% [25]. 


In previous studies of cosmic ray muon events at LEP, QGSJET 01 was used to model hadronic in¬ 
teractions. Apart from the way in which nonlinear effects are modelled, another significant difference 
between this earlier version of the model and QGSJET 11-03/04 is the deeper shower maximum, X max , 
used in the later versions. This results in a steeper lateral muon distribution and an associated increase 
of the muon density close to the core of the shower, which can also have an impact on the observed rate 
of HMM events. 


When generating cosmic ray events, the core of each shower was scattered randomly at ground level 
over an area covering 205 x 205 nr centred upon the nominal LHC beam crossing point. This area was 
chosen to minimise the number of events to be generated without creating any bias on the final results. 
We found that, when the core was located outside this area, only a very small number of events gave 
rise to muons crossing the TPC and these events were always of low multiplicity (/V u < 4). Therefore, 
neglecting these events does not affect the results reported in this paper. 

To have a fast and flexible method of estimating several important parameters and observables involved 
in the analysis, we started with a simplified Monte Carlo simulation. This simulation did not explicitly 
model interactions in the rock above the experiment. Instead, the trajectories of the muons arriving at the 
surface were simply extrapolated as straight lines to the depth of ALICE while imposing an energy cut 
E^> 16 GeV/ cos(0), where 6 is the zenith angle of the muon. All muons passing this cut and crossing 
an area of 17 nr 2 , corresponding to the horizontal cross-sectional area of the TPC, were considered to be 
detected. 


To understand the complete sample of the recorded data, including the origin of low muon multiplicity 
events, we generated events initiated by the interaction of proton and iron ( 56 Fe) primaries with energies 
E > 10 12 eV. This revealed that most single muon events stem from primaries in the energy range 10 12 < 
E < 10 13 eV, while primaries in the energy range 10 13 < E < 10 14 eV produce muon multiplicities 
typically in the range from 1 to 4, independent of the mass of the primary cosmic rays. Primaries 
with energies below 10 14 eV therefore produce a negligible contribution to multi-muon events (/V u > 4) 
that are of interest in this study. Consequently, only energies E > 10 14 eV were considered in the full 
simulation. 

The first step in the analysis was to attempt to reproduce the measured MMD in the low-intermediate 
range of multiplicity (7 < N, < 70). Samples of proton and iron primary cosmic rays were generated 
in the energy range 10 14 < E < 10 18 eV and with zenith angles in the interval 0° < 0 < 50°. The 
composition of cosmic rays in this energy range is a mixture of many species of nuclei in a ratio that 
is not well-known and which varies with energy. To simplify the analysis and interpretation of the data 
we have modelled the primary cosmic ray flux using a pure proton sample, representing a composition 
dominated by light nuclei, and a pure iron sample, representing a composition dominated by heavy 
nuclei. In relation to the MMD, the proton sample provides a lower limit on the number of events for a 
given multiplicity, while the iron sample provides an upper limit. A typical power law energy spectrum, 
E has been adopted with a spectral index y = 2.7 ±0.03 for energies below the knee (Z± = 3 x 10 15 
eV) and y/ ; = 3.0 ±0.03 for energies above the knee. The total (all particle) flux of cosmic rays has 
been calculated by summing the individual fluxes of the main chemical elements at 1 TeV l26ll where 
measurements are most precise. The flux was estimated to be F(1 TeV) = 0.225 ±0.005 (nr 2 s sr TeV) -1 . 

All events generated with energies E > 10 14 eV were subsequently considered for a complete analysis 
using a detailed simulation taking into account all possible interactions in matter surrounding the exper¬ 
iment. In each event, all muons were extrapolated to the horizontal mid-plane of the experiment and 
flagged if they hit an enlarged area of 36 nr 2 centred upon the TPC with no restriction on the energy of 
the muons. All flagged muons were recorded along with their position and momentum at ground level 
and used as input to the ALICE simulation framework. In this framework, the ALICE experimental hall 
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and the environment above and around the apparatus as well as all the detectors are accurately described. 
Flagged muons were propagated through this environment with GEANT3 12711 . Any muon that crossed 
the detector apparatus was treated by a detector response simulation that produced pseudo-raw data, 
which was then processed with the same reconstruction code that was applied to real data, including the 
TPC tracking algorithm and the Lack matching algorithm developed for this analysis. 


4.1 The muon multiplicity distribution 


We generated simulated events equivalent to 30.8 days live time to permit direct comparison with the 
data without the need to apply an arbitrary normalisation factor. A comparison of the trigger corrected, 
measured MMD with the simulations is shown in Fig. [6] For ease of comparison, the points obtained 
with the simulations were fitted with a power-law function to obtain the curves for proton and iron. 


At lower multiplicities, corresponding to lower primary energies, we find that the data approach the 
proton curve, which represents a light ion composition of the primary cosmic ray flux, while higher 
multiplicity data lie closer to the iron curve, representing a heavier composition. The limited statistics in 
the range > 30 does not allow for a precise, quantitative study of the composition but suggests that 
the average mass of the primary cosmic ray flux increases with increasing energy, a finding consistent 
with several previous experiments 281 43 lh . 


The errors in Fig. [6] arc shown separately (statistical and systematic) for data, while for Monte Carlo 
they are the quadrature sum of the statistical and systematic uncertainties. The systematic errors in 
the simulations take into account uncertainties in the flux of cosmic rays at 1 TeV, the slope of the 
energy spectrum below and above the knee, the description of the rock above the experiment and the 
uncertainty in the the number of days of data taking (detector live time). The largest contribution to 
the systematic error is due to the uncertainty in the spectral index below the knee (y = 2.7 ± 0.03), 
which results in an uncertainty of approximately 15% in the MMD. The error in the description of the 
rock above the experiment corresponds to an uncertainty in the energy threshold of the muons reaching 
the detector, which results in a systematic error of approximately 4%. Each of the other uncertainties 
gives a contribution of around 2% to the systematic error. For muon multiplicities A a > 30, statistical 
uncertainties are dominant. 


Following success in describing the magnitude and shape of the MMD over this intermediate range of 
multiplicities (7 < A4 < 70) we have used the same simulation framework to study the frequency of 
HMM events. Since these are particularly rare events, a very high statistics sample of simulated HMM 
events was required to permit a meaningful quantitative comparison. 


4.2 High muon multiplicity events 

Taking the dataset as a whole, corresponding to 30.8 days and a mixture of running conditions, we find 5 
HMM events with muon multiplicities > 100 (as can be seen in Fig. [5J giving a rate of 1.9 x 10 6 Hz. 
Each of these events were examined closely to exclude the possibility of “interaction” events. The highest 
multiplicity event reconstructed in the TPC was found to contain 276 muons, which corresponds to a 
muon areal density of 18.1 m 2 . For illustration, a display of this event is shown in Fig. |7] The zenithal 
and azimuthal angular distributions of the muons from the same HMM event are shown in Fig. [8j while 
the spatial distribution of matched and single-track muons at the TPC mid plane is shown in Fig. [9] We 
note that the majority of single-track muons are reconstructed near the ends of the TPC where muons 
may enter or leave the active volume without producing a track either the upper or lower halves of the 
detector. 

One of the aims with this study is to compare the rate of HMM events obtained from simulations to 
the measured rate. To limit the effect of fluctuations in the number of simulated HMM events, we have 
simulated a live time equivalent to one year with CORSIKA 6990 using QGSJET11-03 for the hadronic 
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Fig. 6: The measured muon multiplicity distribution compared with the values and fits obtained from CORSIKA 
simulations with proton and iron primary cosmic rays for 30.8 days of data taking. The errors are shown separately 
(statistical and systematic) for data, while for Monte Carlo they are the quadrature sum of the statistical and 
systematic uncertainties. 



Fig. 7: Event display of a multi-muon event with 276 reconstructed muons crossing the TPC. 
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Fig. 8: Zenithal and azimuthal distribution of the multi-muon event with 276 reconstructed muons. 
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Fig. 9: Spatial distribution of the 276 recostructed muons indicating matched and single-track muons. 

interaction model. The simplified Monte Carlo used as a first step of the analysis demonstrated that only 
primaries with energy E > 10 16 eV contribute to these events. Therefore, only events in the range of 
primary energy 10 16 < E < 10 18 eV have been generated to achieve an equivalent of 365 days exposure 
for both proton and iron primaries. 

The estimated maximum fiducial area of the TPC due to its horizontal cylindrical geometry and cut 
on the minimum number of TPC space points is 17 ±0.5 nr 2 . The estimated error in the number of 
reconstructed muons, A±, counting both matched and single-track muons, is around 5% for A u > 100. 
HMM events are therefore events with a muon areal density p Ll > 5.9 ±0.4 m 2 and correspond to a 
rate of 1.9 x 10 6 Hz at the underground location of ALICE. Based upon the number of observed HMM 
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CORSIKA 6990 
QGSJET11-03 
Simple MC Full MC 


CORSIKA 7350 
QGSJET 11-04 
Simple MC Full MC 


Run 


proton iron proton 


iron 


proton iron proton 


iron 


40 

61 

27 

51 

41 

72 

30 

52 

40 

64 

24 

42 

42 

88 

32 

71 

31 

43 

25 

31 

48 

78 

29 

62 

26 

52 

20 

34 

46 

84 

35 

61 

33 

64 

22 

53 

36 

83 

31 

58 


Table 1: Number of HMM events for each run obtained with the simplified Monte Carlo and the full simulation. 
Each run is equivalent to 365 days of data taking. The events have been generated using CORSIKA 6990 with 
QGSJET 11-03 and CORSIKA 7350 with QGSJET 11-04. 



CORSIKA 6990 

CORSIKA 7350 


QGSJET 11-03 

QGSJET 11-04 


proton iron 

proton iron 

<N> 

23.6 42.2 

31.4 60.8 

<7 

1.3 (5.5%) 5.0(12%) 

1.1 (3.7%) 3.5 (5.7%) 


Table 2: Mean value and statistical uncertainty in the number of HMM events for 365 days live time calculated 
using the full simulation. 

events, the estimated relative statistical uncertainty is 45%, giving an error in the rate of ±0.9 x 10 6 Hz. 

The rate of HMM events obtained with the Monte Carlo can be compared with the observed rate. Since 
we have simulated samples of HMM events corresponding to one year live time, the statistical uncertainty 
in the simulated rate will be lower than that in the measured rate. Results obtained for the number of 
HMM events expected in one year from both the simplified Monte Carlo and the full simulation (the first 
of five statistically independent simulations) are shown in the first row of Table Q] Comparison of the 
results demonstrates that the detailed modelling of the underground environment has about a 30% effect 
on the number of HMM events. Due to the small numbers of HMM events we reused the same simulated 
EAS sample to perform four additional simulations by randomly assigning the core of each shower over 
the usual surface level area of 205 x 205 nr 2 . Given that the acceptance of the TPC is almost 3000 times 
smaller, this ensures that the samples are statistically independent. A summary of the results obtained for 
all five simulations is presented in Table |T|for both CORSIKA 6990 with QGSJET 11-03 and CORSIKA 
7350 with QGSJET 11-04. 

Final values for the HMM event rate for proton and iron primaries were calculated by taking the average 
value obtained from the five simulations, while the statistical uncertainty was estimated from the standard 
deviation of the 5 values from the mean. Table [2]summarises the mean number of HMM events expected 
in one year for each primary ion calculated from the full simulation. 

There are two major contributions to the systematic uncertainty on the number of HMM events. The first 
contribution stems from the muon reconstruction algorithm. To estimate its contribution we took the first 
simulated sample, corresponding to 365 days of data taking, for each element and each CORSIKA code 
version and redetermined the number of HMM events using different tunes of the track selection and 
matching algorithms. The second contribution stems from the uncertainties of the parameters used in the 
simulations, as discussed in section 14.11 This was estimated to give an uncertainty in the predicted rate 
of HMM events of approximately 20%. Due to the large sample used in the simulations (365 days), the 
systematic uncertainty is dominant, while in the data (30.8 days) the statistical uncertainty is dominant. 
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HMM events 

CORSIKA 6990 
QGSJET 11-03 
proton iron 

CORSIKA 7350 
QGSJET 11-04 
proton iron 

Data 

Period [days per event] 

15.5 

8.6 

11.6 

6.0 

6.2 

Rate [xl(T 6 Hz] 

0.8 

1.3 

1.0 

1.9 

1.9 

Uncertainty (%) (syst + stat) 

25 

25 

22 

28 

49 


Table 3: Comparison of the HMM event rate obtained with the full simulation and from measurement. 

The systematic uncertainties have been added in quadrature to the statistical uncertainty in the final 
comparison of the observed rate of HMM events with that obtained from the Monte Carlo simulations. 


5 Results 


In Table [3] we present the results of this analysis where we compare the rate of simulated HMM events 
with the measured rate. We note that the pure iron sample simulated with CORSIKA 7350 and QGSJET 
11-04 produces a HMM event rate in close agreement with the measured value. The equivalent rate ob¬ 
tained with CORSIKA 6990 and QGSJET 11-03 is lower, although still consistent with the measured rate. 
The difference between the two simulations comes primarily from the hadronic model used to generate 
the EAS. It is more difficult to reconcile the measured rate of HMM events with the simulated rate ob¬ 
tained using proton primaries, independent of the version of the model. However, the large uncertainty in 
the measured rate prevents us from drawing a firm conclusion about the origin of these events, although 
heavy nuclei appear to be the most likely candidates. Therefore, an explanation of HMM events in terms 
of a heavy primary cosmic ray composition at high energy and EAS described by conventional hadronic 
mechanisms appears to be compatible with our observations. This is consistent with the fact that they 
stem from primaries with energies E > 10 16 eV, where recent measurements [32, 33] suggest that the 
composition of the primary cosmic ray spectrum is dominated by heavier elements. 


Finally, we have investigated the distribution of simulated EAS core positions at the location of ALICE 
for each of the HMM events simulated with iron primaries using CORSIKA 7350 and QGSJET II- 
04 in Table Q] equivalent to 5 years of data taking. The distribution is shown in Fig. |T0j where the 
colour of each point indicates the energy associated with the primary cosmic ray so as to give a visual 
representation of the correlation between the distance of the core from the centre of ALICE at surface 
level and the energy of the primary cosmic ray. We note that the shower cores of all HMM events fall 
within an area of approximately 140 x 140 nr 2 centred upon ALICE, which is located at the origin in 
Fig.lTOl The average distance of the shower core from the centre of ALICE for all events is 19 m and the 
RMS value of the distribution is 16 nr. Primaries with an energy E > 3 x 10 17 eV, corresponding to the 
highest energy interval studied in this analysis, produce larger showers that may give rise to HMM events 
when the shower core falls farther from the location of ALICE. In this case, the mean of the shower core 
distribution from the centre of ALICE is 37 nr and the RMS value of the distribution is 18 m. 


6 Summary 

In the period 2010 to 2013, ALICE acquired 30.8 days of dedicated cosmic ray data recording approx¬ 
imately 22.6 million events containing at least one reconstructed muon. Comparison of the measured 
muon multiplicity distribution with an equivalent sample of Monte Carlo events suggests a mixed-ion 
primary cosmic ray composition with an average mass that increases with energy. This observation is 
in agreement with most experiments working in the energy range of the knee. Following the successful 
description of the magnitude of the MMD in the low-to-intermediate range of muon multiplicities we 
used the same simulation framework to study the frequency of HMM events. 
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Fig. 10: The surface level spatial distribution of the cores of simulated EAS giving rise to more than 100 muons in 
the ALICE Time Projection Chamber. The simulation was for iron primaries in the energy range 10 16 — 10 18 eV 
and corresponds to the equivalent of 5 years of data taking. 


High muon multiplicity events were observed in the past by experiments at LEP but without satisfactory 
explanation. Similar high multiplicity events have been observed in this study with ALICE. Over the 
30.8 days of data taking reported in this paper, 5 events with more than 100 muons and zenith angles less 
than 50° have been recorded. We have found that the observed rate of HMM events is consistent with 
the rate predicted by CORSIKA 7350 using QGSJET11-04 to model the development of the resulting air 
shower, assuming a pure iron composition for the primary cosmic rays. Only primary cosmic rays with 
an energy E > 10 16 eV were found to give rise to HMM events. This observation is compatible with a 
knee in the cosmic ray energy distribution around 3 x 10 15 eV due to the light component followed by a 
spectral steepening, the onset of which depends on the atomic number (Z) of the primary. 

The expected rate of HMM events is sensitive to assumptions made about the dominant hadronic produc¬ 
tion mechanisms in air shower development. The latest version of QGSJET differs from earlier versions 
in its treatment of forward neutral meson production resulting in a higher muon yield and has been re¬ 
tuned taking into account early LHC results on hadron production in 7 TeV proton-proton collisions. 
This is the first time that the rate of HMM events, observed at the relatively shallow depth of ALICE, 
has been satisfactorily reproduced using a conventional hadronic model for the description of extensive 
air showers; an observation that places significant constraints on alternative, more exotic, production 
mechanisms. 


Compared to the previous studies at LEP, there are two distinguishing aspects of this work that have led to 
these new insights into the origin of HMM events. The first has been the ability to generate large samples 
of very energetic cosmic rays, allowing for a more reliable estimate of the expected rate of these events. 
The second, and more important, aspect has been the recent advances in the hadronic description of EAS. 
This is a continually evolving field. We note that in a preparatory study [18] earned out by ALICE in 
2004, using an older version of CORSIKA (version 6031), no HMM events were observed in the MMD 
distribution simulated for 30 days of data taking with a pure iron primary cosmic ray composition. In the 
present work, Table [3]gives a quantitative comparison of the rate of HMM events predicted by two more 
recent versions of CORSIKA and QGSJET, illustrating the evolution of the hadronic description of EAS 
in recent years. Only in the latest version of the model there has been a significant increase in the rate of 
HMM events that better approaches the rate observed in this study. 
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